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Abstract 

We consider a variational problem for three-dimensional (3D) classical lattice models. 
We construct the trial state as a two-dimensional product of local variational weights that 
contain auxiliary variables. We propose a stable numerical algorithm for the maximization 
of the variational partition function per layer. The numerical stability and efficiency of the 
new method are examined through its application to the 3D Ising model. 

1 Introduction 

The density matrix renormalization group (DMRG) has been applied to a variety of one- 
dimensional (ID) quantum systems and 2D classical systems PJIBIE]- This method has also 
been applied to finite-size 2D quantum systems, which can be represented as ID quantum 
systems with long-range interactions. [3J [2] Despite this success, no extension of the DMRG 
to infinitely large higher-dimensional systems has been reported to this time. This is partially 
because the density matrix eigenvalues decay very slowly in higher-dimensional systems, and the 
RG transformation, which is obtained from the diagonalization of the density matrix, becomes 
ineffective. (S] 

The numerical efficiency of the DMRG for ID quantum and 2D classical systems is due 
in part to its variational background, [U Hj in which the trial state is constructed as a 
product of orthogonal matrices that represent the block-spin (or the renormalization group) 
transformations. Such a construction of the variational (or trial) state can be generalized to 
higher dimensions. A simple example is to use a 2D classical system as a variational state for 
2D quantum and 3D classical systems. Martm-Delgado et al. employed the 6-vertex model as 
a trial wave function for 2D lattice spin/electron systems. :9] Okunishi and Nishino considered 
the direct extension of the Kramers- Wannier approximation jlUj to the 3D Ising model, in which 
the variational state is the 2D Ising model subject to an effective magnetic field. These 
examples have demonstrated the potential usefulness of the 2D variational state, constructed 
as a product of local variational weights. 

For the purpose of obtaining a better variational result, Nishino et al. developed a numerical 
method, the tensor product variational approach (TPVA), which automatically optimizes the 
2D variational state using the solution of a self-consistent equation. ^2] They applied this 
method to both the 3-state (q = 3) Potts and the Ising models on simple cubic lattices, and 
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Figure 1: The Boltzmann weight Wjj {o"|cr} and the transfer matrix T[<t|(t] in the case N = 2. 

reported that TPVA yields a better estimate of the transition temperature for the q = 3 Potts 
model. This is because the number of variational parameters in the trial state for the Potts 
model, 3 4 , is greater than that for the Ising model, 2 4 . 

In general, the TPVA gives a lower variational free energy when the trial state contains a 
larger number of parameters. A way of increasing this number is to employ a variational state 
that contains auxiliary variables, for example block spin variables. ^SI 1 This generalization, 
however, introduces a serious instability into the numerical optimization of the variational 
state. [T3] 

In this paper we report on a method of stabilizing the numerical optimization process. We 
introduce an orthogonal condition placed on a small change of the local variational weight. 
This condition prevents an 'unexpected decrease' of the norm of the variational state, which 
caused numerical instability in previous calculations. ^3] I n the next section, we review the 
formulation of the TPVA. jl3U12llT3] In §3 we consider the stability of the optimization process. 
The numerical efficiency of the stabilized TPVA is examined in §4 in the case that it is applied 
to the 3D Ising model. Conclusions are given in §5. 



2 Tensor product variational approach 

We consider the 3D Ising model on a simple cubic lattice as an example. The system size is 
2N x 2N x oo in the X-, Y-, and Z-directions. On each lattice point there is an Ising spin a, 
with a = dbl. A ferromagnetic interaction —Jaa' exists between nearest-neighbor spins. 

As shown in Fig. ^ the transfer matrix T connects adjacent spin layers [a] and [a], where 
[a] represents all the spins in a layer: 

[a] = {a xl , a 2 i, • • • , (J 2 N 1> °12> °2 2i ' " " > &2N 2N J ■ (1) 

For simplicity, we consider a symmetric transfer matrix, 

2JV-1 

T[a\a]= J] W% j) {a\a}, (2) 

where Wg J ' {ajcr} represents a local Boltzmann weight with respect to a unit cube at the 
position (A, Y) = (see Fig. ^). We denote a spin plaquett, a group of the nearest 4 spins, 
as 

W} = ( a ij a i+ij a i+ij+l) = ( a ij a ij> a i'j ( 3 ) 

There is another way of tuning the local variational weight by way of the density matrix renormalization in 
vertical direction (see Ref. |14p. 
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Figure 2: A graphical representation of the local tensor and the construction of the trial 
state VP in the case N = 2. The circles and ovals denote the Ising spins a and the n-state 
auxiliary variables £, respectively. 



using the i' = i + 1 and f = j + 1. With this notation, we can write the local Boltzmann 
weight of the 3D Ising model as follows: 



W£ j) {a\a} = exp 



J 



4k B T 



<T ij a i'j + a i'j a i'j' + a i'j' a ij' + a ij' a ij + a ij a i'j 



+ a i'j (J i'j' + Ci'j'Vi? + a ij' a ij + a ij a ij + a i'j a i'j + a i'j' a i'j' + a ij ,(T ij' 

We have thus expressed the 3D Ising model as a special case of the 'interaction round a face' 
(IRF) model. 

For an arbitrary variational function the Rayleigh ratio 

E <S>[a}T[a\a]*[a} 
Xm = m s ffl (4) 

gives the variational partition function per layer. In the framework of the TPVA, the trial state 
^ is constructed as a contracted product of local variational weights as 

* m = e n W y ? j y °i y - e n ^ ( : ) - 

where the variables ^ ■ of the local weight V are auxiliary variables that can be in one of n 
states (0,1, ... ,n— 1), while the spin variables ■ can be in one of 2 states (±1) 2 (see Fig.0. 
We have expressed a group of auxiliary variables as [£], analogously to [cr] in Eq. Q. 3 We arc 
interested in the bulk properties of lattice models, and therefore we consider the case in which 
the system size 2N is sufficiently large and assume that the local variational weig ht v( y ) is 
position independent. Hereafter, we refer to this variational state as the 'tensor product state' 
(TPS). 

Because both the trial state $ and the transfer matrix T are written as products of local 
variational weights, both the numerator and denominator of Eq. @ can also be expressed as 
products of stacked local weights: 

<*m*> = E llV^ (III) W^{a\a] VM (^j) , 

[elite] 1 M 



2 Two kinds of tensor products are known. One is the IRF type |12j and the other is the vertex type. We 
use the former throughout this section. 

3 When n = 1, the variational function $ does not contain any auxiliary variables. 
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(I'll') 



(6) 



In other words, both (\I'|\I') and (\E , |7~|\I r ) are partition functions of stacked two-dimensional 
classical systems. This structure enables us to calculate and (^IT^) accurately by 

use of numerical renormalization techniques. |16| I17 [ IH] Thus, the variational partition function 
A^) can be easily calculated for an arbitrary TPS constructed from the local variational weight 
V. 

We maximize A(\I/) by tuning elements of the local weight V. In order to clarify the 
variational structure with respect to V, let us divide (^l^} into two parts: (i) the adjacent 
local weights V^ NN ^ and V^ NN ^ at the center, and (ii) the rest, consisting of the stacked classical 
system with a puncture at the center ^Sj, expressed as 

>} {a } )-S n v<« (%)vm (W 



A 



[f YA& 



(7) 



where the configuration sums are taken over all values of the variables a and £, except for those 



(i,j)^(N,N) 



indicates that 



at the center that are variables of V( Nir > and V^ NN \ The notation n 

y( NN ) i s no t included in the product. In the same manner, we divide (*|T|*> into (i') V^ NN ^> 
and \f( NN \ and (iii) the rest, expressed by 



B 



{*} 
{0 



M 



w, 



(NN) 



{a\a} J2 II y{ij) 

r]> (i,j)^(N,N) 



wg j) {a\a}V^ 



M 



(8) 



which corresponds to a partially punctured stacked classical system. The meanings of the 
matrices A and B are represented graphically in Fig. |21 Using this new notation, we can 
rewrite the Rayleigh ratio, Eq. (jlj, as the following ratio of bilinear forms: 



A(*) 



„E/™»( { «-i) B Gf 

{£},{(} 



} 

{£} 



y(NN) 



{£}) 



£ V^(\f } )A(\f } \U)vi^(^) -WAV)' 

{?},{?} 



(9) 



Here, we have interpreted the variational weig ht at the center y (AW) as a (2n) 4 -dimensional 
vector and have written it simply as \V). 
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As we have assumed that the system size 2N is sufficiently large, the variation of A(^) with 
respect to a uniform modification of local weights, 

tfA(¥) = ]^7) 6v{ij) + [ higher order terms l ' ( 10 ) 

ij 

is almost proportional to dX(tp)/dV^ NN \ which is the contribution from the local variation 
only at the center. |121 I13| From the optimal condition dX(i^)/dV^ NN ^ = 0, we obtain the 
generalized eigenvalue problem 

B\V) = XA\V). (11) 

We use this equation when we optimize the TPS. Note that Eq. ()11|) is a non-linear equation 
with respect to the local weights V, since the (2n) -dimensional matrices A and B themselves 
are functionals of V^ NN \ Therefore, Eq. (jlljl should be solved self-consistently by use of 
successive, small improvements of the local variational weig ht vw. 

In Ref. |13j . Nishino et al. apply the TPVA to the 3D Ising model, which is represented 
as a symmetric 16-vertex model. For the case n = 2, they succeeded in optimizing the TPS. 
However, they encountered numerical instability in the case n > 3. The reason for this is that 
the matrix A is not always positive definite during the process in which the variational state 
is being improved, although A should be positive definite after the optimization is completed. 
We discuss the cause of this instability and propose a method of stabilization in the following. 



3 Stabilization of the self-consistent improvement process 

Consider an infinitesimal change of the local weight at the center of the system, 

\V) - \V) = \V) + e\X), (12) 
where \X) is an arbitrary (2n) -dimensional vector. The Rayleigh ratio, Eq. Q, is modified as 



A' 



(V'\B\V) _ (V\B\V) 
(V'\A\V') ~ (V\A\V) 



( (Y\B\X) r2 (X\B\X) \ 

(V\B\V) (V\B\V) 
(V\A\X) j(X\A\X) 



(13) 



(V\A\V) (V\A\V)J 



If e is sufficiently small, we can expand A' as 

1 + 2s((3 - a) + 0(e 2 ] 

with 



A 



A' = A 
(V\B\V) 



a 



(V\A\X) 



(v\b\x) 



(14) 



(15) 



(V\A\V) ' (^1^1^) ' r (V\B\V) 

It seems appropriate to select \X) to maximize — a. However, such a choice of \X) tends to 
cause the expectation value of the denominator (V'|«4|V) to decrease, and after several self- 
consistent iterations, (V|«4|V) becomes very small or negative. This is a cause of numerical 
instability in the previous TPVA algorithm. We prevent this anomalous decrease of (V|^4|V) 
by choosing \X) to satisfy 

(V'\A\V) = (V\A\V) + 0{e 2 ) . (16) 
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This is equivalent to choosing \X) to be orthogonal to .A|V), which gives a = 0. With such a 
choice, the maximization of A is carried out with respect to (5 only. A reasonable choice of \X) 
can be obtained using the Schmidt orthogonalization: 

\X)=B\V)-A\V){V\AB\V). (17) 

The \X) thus obtained at least yields non- negative (3. We include \X) in the self-consistent 
calculation of the TPVA as follows: 

(1) Prepare an arbitrary initial variational weight |V). 

(2) Calculate the matrices A and B using a numerical RG method. 

(3) Obtain \X) from Eq. (fT7|) . 

(4) Improve the local weight according to \V) = \V) + s\X). 

(5) Terminate the modification of |V) if the computed thermodynamic functions have con- 
verged; otherwise go back to step (2). 

The third step is the main difference between the present and the previous TPVA algorithms. 
Since (V|.A|V) changes only by an amount of order e 2 at each iteration, a large number of 
iterations is necessary to realize convergence, especially when e is small. 

4 Numerical results 

We now confirm the numerical stability of the new algorithm in the case that it is applied to 
the 3D Ising model. Hereafter, we set J/k B = —1 and denote by m the number of the states 
of the block-spin variable used in the numerical RG calculation |lfi[ I17j , which is used for the 
calculation of A and B. 

Figure 0] shows the convergence of the spontaneous magnetization (a) with respect to the 
number of numerical iterations at the temperature T = 4.504. (The selected temperature 
T = 4.504 is slightly lower than the critical temperature T c = 4.5115 obtained by Monte Carlo 
simulations ^S] which is considered to be the most reliable result.) The parameter e in Eq. 1)12(1 
is chosen to be 10 -3 for normalized \V) and \X). In both the cases n = 1 and n = 2, (a) begins 
to converge monotonically to the final value after several hundred iterations. 

Figure Efa) displays (a) thus calculated for n = 1 (no auxiliary variables) and n = 2 with 
m = 5. In the region T < 4.2, the difference between these cases is not visible. As shown in the 
inset, near the critical temperature, the calculated (a) with n = 2 decays more rapidly than 
that with n = 1. The estimated transition temperatures, where the obtained (a) falls to zero, 
are T = 4.57 (n = 1) and 4.55 (n = 2). 

To this point, we have expressed the 3D Ising model as a 3D IRF model. The Ising model 
can also be expressed as a special case of the symmetric 64-vertex model. ^Hl Applying the 
stabilized TPVA algorithm to this vertex expression, we obtain the (a) plotted in Fig. E^b). 
In this case, the estimated transition temperatures are T = 4.533 (with n = 2 and m = 16) 
and 4.525 (with n = 3 and m = 12). All of these calculated transition temperatures are higher 
than T C MC = 4.5115, obtained by Monte Carlo simulation, R is thus seen that with the 
TPVA, the ordered phase tends to be stabilized. 

We finally compare the stabilized numerical algorithm considered here with the previous 
one. |13J The stabilization provided by the Schmidt orthogonalization, expressed by Eq. (|17|). 
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Figure 4: Convergence of the spontaneous magnetization (a) with respect to the number of 
numerical steps at T = 4.504. 

enables us to perform the calculation for those cases n = 2 (IRF expression) and n = 3 
(Vertex expression) for which the previous algorithm does not yield a convergent result. With 
regard to the computational time required to obtain the convergent numerical result for the 
cases n = 1 (IRF expression) and n = 2 (Vertex expression), for which the previous and the 
stabilized algorithms give the same numerical result, the stabilized algorithm is about 10 times 
slower than the previous algorithm. In the variational calculation employing the TPVA, the 
numerator (V|£?|V) should be maximized and the denominator (V|«4|V) should be minimized 
through successive improvements of V. Though (V|2?|V) increases rather rapidly, (V|.4.[V) 
decreases slowly, because the stabilization condition Eq. Q16JI restricts the size of the change of 
C^|-4|V) to order e 2 . 

5 Conclusion and discussion 

We have proposed a stabilized partition function maximization process for numerical calcula- 
tions employing the TPVA, imposing the orthogonality condition on the small change made to 
the local variational weight. This improvement enables us to obtain the magnetization using a 
TPS that contains n-state auxiliary variables, in particular, with n = 2 for an IRF-type TPS 
and n = 3 for a vertex- type TPS. 

The orthogonality condition expressed by Eq. (|17|) stabilizes the numerical calculation, but 
it also slows the convergence to the variational maximum, because it causes the improvement 
of the denominator of the variational formulation to be very slow. Contrastingly, the vertical 
density matrix approach (VDMA) ^1], which also uses a TPS as its trial state, exhibits rapid 
convergence. In the VDMA, the local weight is improved through a RG transformation ob- 
tained from the diagonalization of the density matrix, and thus both the denominator and the 
numerator of the variational ratio are improved equally rapidly. The VDMA, however, has the 
shortcoming that it does not fully improve the TPS, because the RG transformation used in 
the VDMA is not determined so as to maximize the variational partition function per a layer. 
In higher dimensions, the direct diagonalization of the density matrix does not give the most 
appropriate RG transformation. Our next project is to combine the advantages of the TPVA 
and the VDMA. 




1 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 ' 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 

T[MJ TIM J 



Figure 5: Calculated spontaneous magnetizations, (a) (a) obtained for (n,m) = (2,5). The 
curve here was drawn using the Spline interpolation. The inset displays (a) near the criticality 
for (n,m) = (1,5) (circles) and (2,5) (triangles), (b) (a) obtained from the TPVA applied to 
the Ising model represented as the 64- vertex model |13| for (n,m) = (2, 16). The inset plots 
the data for (n,m) = (2, 16) (circles) and (n, m) = (3, 12) (triangles). 
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